function [Ees,Eep,fi1,fi2,fi3]=f_2mlyrAmp(n0,n1,k1,n2,k2,n3,d1,d2,m,fi0,lambda)

%calculate transmittance and reflectance of a multi-layer system

%Written by Feiran Huang 05/31/99. Last Update by JW Haus
%              08/13/2002 Amplitudes of the fields added.

%-------------------------

% 





   fi1=asin(sin(fi0)*n0/(n1-i*k1));
phi1=fi1;
   fi2=asin(sin(fi0)*n0/(n2-i*k2));
phi2=fi2;
   fi3=asin(sin(fi0)*n0/n3);
phi3=fi3;


   [r01p,t01p,r01s,t01s]=f_rtamp(n0,0,n1,k1,fi0,fi1);
   [r12p,t12p,r12s,t12s]=f_rtamp(n1,k1,n2,k2,fi1,fi2);
   [r21p,t21p,r21s,t21s]=f_rtamp(n2,k2,n1,k1,fi2,fi1);
   [r23p,t23p,r23s,t23s]=f_rtamp(n2,k2,n3,0,fi2,fi3);

   delta01=0;

   delta12=2*pi*(n1-i*k1)*d1*cos(fi1)/lambda;
   delta21=2*pi*(n2-i*k2)*d2*cos(fi2)/lambda;
   delta23=2*pi*(n2-i*k2)*d2*cos(fi2)/lambda;

   C01p=[1 r01p; r01p 1];

   C12p=[exp(i*delta12) r12p*exp(i*delta12); r12p*exp(-i*delta12) exp(-i*delta12)];

   C21p=[exp(i*delta21) r21p*exp(i*delta21); r21p*exp(-i*delta21) exp(-i*delta21)];

   C23p=[exp(i*delta23) r23p*exp(i*delta23); r23p*exp(-i*delta23) exp(-i*delta23)];

   

   
   %   field amplitude calculations p-polarization
   
   % transmitted field amplitude
   Temp=[1 0]';
   Eep(:,2*m+2)=Temp;
   % Amplitude in the last layer
   Temp=C23p*Temp/t23p;
   Eep(:,2*m+1)=Temp;
   for j=0:m-2
      Temp=C12p*Temp/t12p;
      Eep(:,2*(m-j))= Temp;
      Temp=C21p*Temp/t21p;
      Eep(:,2*(m-j)-1)= Temp;      
   end
   Temp=C12p*Temp/t12p;
   Eep(:,2)=Temp;
   Temp=C01p*Temp/t01p;
   Eep(:,1)=Temp;
   
   % NORMALIZING AMPLITUDE  
    Eep=Eep/Eep(1,1);


   
   C01s=[1 r01s; r01s 1];

   C12s=[exp(i*delta12) r12s*exp(i*delta12); r12s*exp(-i*delta12) exp(-i*delta12)];

   C21s=[exp(i*delta21) r21s*exp(i*delta21); r21s*exp(-i*delta21) exp(-i*delta21)];

   C23s=[exp(i*delta23) r23s*exp(i*delta23); r23s*exp(-i*delta23) exp(-i*delta23)];

   
     
 %   field amplitude calculations s-polarization
   
   % transmitted field amplitude
   Temp=[1 0]';
   Ees(:,2*m+2)=Temp;
   % Amplitude in the last layer
   Temp=C23s*Temp/t23s;
   Ees(:,2*m+1)=Temp;
   for j=0:m-2
      Temp=C12s*Temp/t12s;
      Ees(:,2*(m-j))= Temp;
      Temp=C21s*Temp/t21s;
      Ees(:,2*(m-j)-1)= Temp;      
   end
   Temp=C12s*Temp/t12s;
   Ees(:,2)=Temp;
  Temp=C01s*Temp/t01s;
   Ees(:,1)=Temp;

% NORMALIZING AMPLITUDE  
    Ees=Ees/Ees(1,1);

















